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Abstract 

The effect of a linear adsorption isotherm on the onset of fingering instability in a miscible 
displacement in the application of liquid chromatography, pollutant contamination in aquifers 
etc. is investigated. Such fingering instability on the solute dynamics arise due to the miscible 
viscus fingering (VF) between the displacing fluid and sample solvent. We use a Fourier pseudo- 
spectral method to solve the initial value problem appeared in the linear stability analysis. The 
present linear stability analysis is of generic type and it captures the early time diffusion dominated 
region which was never expressible through the quasi-steady state analysis (QSSA). In addition, 
it measures the onset of instability more accurately than the QSSA methods. It is shown that 
the onset time depends non-monotonically on the retention parameter of the solute adsorption. 
This qualitative influence of the retention parameter on the onset of instability resemblances with 
the results obtained from direct numerical simulations of the nonlinear equations. Moreover, the 
present linear stability method helps for an appropriate characterisation of the linear and the 
nonlinear regimes of miscible VF instability and also can be useful for the fluid flow problems with 
the unsteady base-state. 
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I. INTRODUCTION 


The displacement process in porous media has enormous importance in the field of fluid 
dynamics. It features viscous fingering (VF) instability when a more viscous fluid is displaced 
by a less viscous one p]. A confined geometry of Helc-Shaw cell [2] is generally used to 
experimentally observe VF in a homogeneous porous medium. Hornsy presented an insightful 
review on VF instability in both miscible and immiscible fluids [I]. Due to its application 
in many industrial and environmental processes, such as the secondary oil recovery from 
porous rocks [3], pollutant contamination in underground aquifers [3], etc., this problem 
has drawn attention of multidisciplinary researchers for many decades. The objectives of 
these researchers are multiple, such as, understanding morphological instability leading to 
interfacial pattern formation [SHE], to present a suitable stability analysis of hydrodynamic 
instabilities with unsteady base-state [SHU, and many more. VF is also observed in liquid 
chromatography, a flow based separation method in which a given fluid (called the displacing 
fluid) displaces a miscible sample consisting of a solvent and a mixture of dissolved solutes 
(called analyte) [15 j. 

Mathematical challenges of performing linear stability analysis (LSA) for this type of 
hydrodynamic instability is that the base state is not translated unchanged along the flow 
direction as time elapses. Diffusion relaxes the interface between the two fluids, resulting 
this as an unsteady base-state problem. The challenge is to capture the linearly unstable 
modes and their onset of instability by suitably incorporating the time evolution of the 
base-state. In order to overcome this difficulty Tan and Hornsy ra used a quasi-steady 
state approximation, which assumes that the growth rate of the disturbances is faster than 
the rate of change of the base state. Following Tan and Hornsy |10j . Rousseaux et al. [16] 
performed a linear stability analysis to capture the onset of instability in liquid chromatog¬ 
raphy column using QSSA. Although QSSA method successfully measures the growth rate 
of the perturbations with a certain degree of accuracy, all the perturbations are found to 
be unconditionally unstable for arbitrary small time t > 0. Therefore, QSSA method fails 
to fulfil the most important goal of stability analysis, which is to appropriately measure 
the onset of instability. Growth rate measured from the QSSA method is valid for large 
time only. Thus, one needs to look for an alternative method for linear stability analysis 
to predict the onset of instability accurately. In this direction Ben et al. [9] performed an 
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LSA using spectral analysis and showed an unconditional stability of the system at early 
times. Recently, Kim |T2], Pramanik and Mishra [13] applied QSSA method in a self-similar 
domain to calculate the onset of VF instability in miscible slices. These authors also found 
unconditional stability at early times, which is in accordance with the spectral analysis [9j. 
Contrary to classical single interface VF [10], the base state concentration profile for the 
displacement of a finite slice is not self-similar. Apart from the QSSA method and spectral 
analysis, another way to perform an LSA of problems with an unsteady base state is to solve 
the linearized equations as an initial value problem (IVP) [TD] • The obstacle in this method 
is that the representative initial condition is not known a priori and often the ‘white noise’ 
or random perturbation is used as the initial condition. The main difficulty arises with such 
an initial condition is that it introduces perturbation to the entire system. One of the aims 
of this paper is to present a generic linear stability analysis for flow with unsteady base-state 
such that it captures all the physics accurately and also is computationally more efficient. 

Adsorption of solute on the porous matrix influences the separation in liquid chromatog¬ 
raphy and pollutant contamination in aquifers mm- A theoretical model incorporating 
linear adsorption of the species ruling the dynamic viscosity of the solution in a miscible dis¬ 
placement was analyzed by Mishra et al. [18]. These authors presented an LSA using QSSA 
method and also direct numerical simulations (DNS) of the nonlinear equations. Their study 
revealed that in the presence of a retention governed by V, the effect of solute is similar 
to that of an unretained solute with the dynamic viscosity reduced by a factor of 1 + V. 
To incorporate the dynamics of the carrier fluid and of both the solute and the solvent, 
Mishra et al. [19] proposed a theoretical model with viscosity being ruled by sample con¬ 
centration and solute being a passive scalar adsorbed linearly. They performed DNS using 
a Fourier pseudo-spectral method to study the effect of VF on spatio-temporal distribution 
of the retained solute and the influence of the retention parameter on the fingering of solute 
concentration. Recently, Mishra et al. m, Rana et al. [211 122] numerically studied the 
influence of solvent modulated adsorption on the propagation dynamics of the solute in the 
absence and presence of viscosity contrast, respectively. 

In this context, we perform a linear stability analysis using a Fourier pseudo-spectral 
method to have a more comprehensive understanding about the influence of solute adsorption 
on the onset of instability. For simplicity, we consider the retention to be independent of the 
solvent concentration. Unlike modal analysis, the present LSA calculates separate growth 
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rates for the perturbations associated with each of the physical quantities. It also captures 
the diffusion dominated region at the early time and distinguishes the linear and nonlinear 
regimes, which were never achieved through the existing LSA methods for problem with the 
unsteady base-state. We show that the linear unstable modes for the sample solvent remain 
unaffected by the adsorption of the solute on the porous matrix. It is further identified that 
the onset of fingering instability of the adsorbed solute has a non-monotonic dependence on 
the retention parameter. These results are consistent with the direct numerical simulations 
of Mishra et al. [HJ]. 

The paper is organized as follows. We present the mathematical model of the problem in 
Sec. [TTJ followed by the stability analysis and numerical method of solution of the problem 
in Sec. |III[ Secs. ITVl and M discuss the obtained results without and with the adsorption 


of the solute concentration on porous matrix, respectively. In these sections, the LSA and 


DNS results are compared, followed by concluding remarks in Sec. VI 
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FIG. 1: Schematic of the flow configuration for three component model with coordinate 

axes. 


II. MATHEMATICAL MODEL 

Consider a uniform rectilinear displacement of a sample solvent of width W and viscosity 
/i 2 injected at initial time t — 0 in a two-dimensional, homogeneous, horizontal porous 
medium (or a Hele-Shaw cell) by a carrier fluid or eluent of viscosity /ij (< fi 2 ) (see Fig. [Tj) . 
The permeability of the porous medium is assumed to be constant K, which is equivalent to 
b 2 /12 for a Hele-Shaw cell consisting of two parallel plates separated by a small gap b <C H, 
where H is the width of the Hele-Shaw plates. The sample consists of a solute or analyte 
of concentration c a = c aj2 dissolved in a solvent of concentration c = c 2 . The sample solvent 
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is different from the carrier fluid in which the sample solvent concentration c = 0 and the 
solute concentration is c a = 0. The eluent is injected with a uniform velocity U along the x 
direction from left to right as shown in Fig. [lj 

Assuming that the fluids are neutrally buoyant, incompressible and the dispersion is 
isotropic we can describe the above-mentioned three component model using the non- 
dimensional equations in a Lagrangian frame of reference moving with the velocity U [ 19] , 


V ■ n = 0, (1) 

Vp = ->n(c)(u + e x ), (2) 

— + n • Vc = V 2 c, (3) 

f)c 

(1 + «') + {u- K'e x ) • Vc a , m = V 2 c a , m , (4) 

where p is the dynamics pressure, u = (u, v) is the two dimensional gap-averaged velocity in 
the Lagrangian frame of reference, fi(c) is the dynamics viscosity of the fluids, and e x is the 
unit vector in the x direction. Further, c and c a>m correspond to the solvent concentration 
and the mobile phase solute concentration, respectively. For the non-dimensionalization we 
use U,D/U, and D/U 2 , respectively, as the characteristic velocity, length and time. Here D 
corresponds to the isotropic dispersion tensor. The reference concentration for the solvent 
and the solute concentration are taken as C 2 and c a ^j (1 + fi/), respectively. The log-mobility 
ratio R = ln(/r 2 /pi), retention parameter k' and the dimensionless sample width w are the 
three dimensionless parameters of the problem. The viscosity and concentration are related 
by an Arrhenius type relationship [if], i.e., /i(c) = e Rc . The initial and boundary conditions 
associated with Eqs. 0-0 are: 


u(x,y,t = 0) = (0,0), 


(5) 

[ 1, 

0 < x < w 


c(x,y,t = 0) = Ca, m {x,y,t = 0) = < 


(6) 

o, 

otherwise, 



and 


n = (0,0), ^ (O’ 0 )’ as M ->■ oo, 

dc dc am dv 

= = °> 
ay ay ay 


( 7 ) 

( 8 ) 
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respectively, where w corresponds to the dimensionless width of the sample. The velocity 
boundary condition in Eq. (j8j) corresponds to constant pressure at the spanwise boundaries, 
where the streamwise velocity component u takes arbitrary value [ 23 j. 


A. Stream function formulation 


For a two dimensional flow, the conservation of mass is satisfied by introducing the stream 
function, p)(x,y,t), such that the velocity components are given by u = p> y , v = —p ) x . Now 
taking the curl of the non-dimensional Darcy’s equation (Eq. ([2])) the pressure is eliminated, 
and we obtain 


v 2 ip 

dc 

777 + 


dc 

-R[^+ Q~ y e y 


Vc, 


dp) dp)\ 2 

V -^>'Vc=V-c. 


. t\dc a 

(l + «) 


dt 


+ 


( dp) , dpj\ 
\dy K ' dx) 


V c a ^ m 


= V 2 C 


a,mi 


(9) 

( 10 ) 


( 11 ) 


where e y is the unit vector in the y direction. The initial condition and the longitudinal and 
transverse boundary conditions corresponding to the velocity are expressed in terms of the 
stream function as, 


dp) dp) 

dy dx 


at t — 0, 


^ kb = (0 0) 

dy' dx 1 [ ’ J ’ 


as |x| —>• oo, 


d 2 P) 

dydx 


= 0, \/x, 


( 12 ) 

(13) 

(14) 


respectively. 


III. STABILITY ANALYSIS AND NUMERICAL SOLUTIONS 

In this section, the unsteady base state of the model is presented, followed by the deriva¬ 
tion of the linearized perturbation equations. Further, the numerical solutions of the linear 
stability problem as well as the direct numerical simulations of the fully nonlinear problem 
using a highly accurate Fourier pseudo-spectral method are presented. The growth rate and 
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the hence onset of instability is obtained by projecting the governing equations as an initial 
value problem. 


A. Base state 


For base-state flow we assume u h = (0,0) which implies ^ to be constant. We also assume 
that the base-state concentration is homogeneous in the y-direction. i.e., c b = c b (x,t ) and 
c a,m. = c a,m( x it)- Using Fourier transform and Eq. (|5]), the base-state flow can be written 
in terms of decaying error function solution of step-like initial concentration profiles for the 
solvent and solute, 


c b (x,t) = ^ erf( 


x 


\2y/t 


erf 


x — w 
2 y/t 


(15) 


cb a,m( X ,t) = ~ 


erf 


—erf 


x + n't/( 1 + k') 

2 \ A /(1 + K ') 

x + n't/(l + k!) — w 

2\A/(1 +7) 


(16) 


respectively. It is clear from Eqs. (15) and (16) that the basic state of the stability problem is 
a diffusing front which is not stationary. To enable modal stability theory, the disturbance 
is often decomposed as Fourier modes in the transverse direction. In this process each 
Fourier mode is also assumed to be decoupled which contradicts the empirical observation 
that the concentration disturbances are localized in the downstream direction within the 
diffusive layer. We know that the concentration gradient at the front scales as \ft for one 
that begins as a step profile. Hence, unless the front has been allowed to diffuse initially 
before displacement begins, the rate of change of the base-state can be arbitrarily large for 
£ 1. In a self-similar coordinate, the base state equation for sample solvent c b 

becomes 






(17) 


where £ = x/y/t is the similarity variable. In contrast to the single interface, the base state 
in the present case is not self-similar i.e., c b is dependent on both £ and t and the transient 
effect of the base state can be no more negligible. Thus to avoid such barriers in modal 
linear stability analysis, QSSA is not invoked in this paper. 
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B. Linearized perturbed equations 


The principle of an LSA is to observe whether infinitesimal perturbations introduced to 


the equilibrium state (Eqs. (15) and (16)) amplify or decay in time. For this, perturbations 
to the base-state flow, such that c(x,y,t ) = c b (x,t ) + 9d(x,y,t), c a: m{x,y,t ) = c b am (x,t) + 
9c' am {x,y,t) and p(x,y,t) — p b + 9p'(x,y,t), etc. are introduced. We already mentioned 
in Sec. Ill A that p b is constant, which can be assumed to be equal to zero without loss 
of generality. On substituting these expressions Eqs. 0-0 can be written in terms of 
perturbation quantities as, 

ax ox dy 


9 2 R 


dd dp' dc' Dili' 
dx dx dy dy 


= 0, 


(18) 


dd dc b dtp' 
dt dx dy 

+ 9 2 


vV 


dc' dip' dd dp>' 


dd 

. /\ a pm 


dx dy dy dx 

J 


= 0, 


(19) 


, dc a,mdp' ,dd , 

H—W—^ k —7T V c, 


9 2 


dx dy 
dc' m dp' 


dx 
dc' dp' 


2J 
a.m 


= o. 


( 20 ) 


dx dy dy dx 

In an LSA the perturbations are required to be infinitesimally small, which can be obtained 
by assuming \9\ <C 1. Neglecting the terms containing 9 2 and equating the coefficients of 9 
equal to zero, we obtain the corresponding linearized equations in terms of the perturbation 
quantities, 


vy + r 


f dc b dp' dc' \ 

\ dx dx ^ dy J 


dc' dc b dp' 
dt dx dy 


V 2 c' = 0, 


dd. 


/ 

a.m 


. r dc a,mdp' dd 
+ 0 — -x-1- A- 


- wv = 0, 


( 21 ) 

( 22 ) 

(23) 


dt dx dy dx 

where 6 = 1/(1 + k') and A = + k'). The boundary conditions associated with the 

equations in perturbation quantities are, 


On C.mdCna, ;</./) = (0,0,0), at x = 0 and L„, 
(c',c ' a . m = (0,0,0), at y = 0 and L„ 


(24) 

(25) 

































where L x and L y represent the dimensionless length and width of the computational domain, 
respectively. In order to understand the instability phenomenon we use different initial 


conditions, and it will be described in Sec. HID 


C. Initial value problem and Fourier pseudo-spectral method 


The traditional approach for studying the LSA for unsteady base state is by frozen-time 
approach which is known as QSSA. In the present case, QSSA approach can not predict 
the growth rate of individual flow variables (see Appendix B 1[ ), and hence it can not meet 
our principal aim i.e., analysing the evolution of growth function associated with solute 
concentration c' a m in relation to the growth function of solvent concentration c !. We present 
a linear stability analysis, which carefully controls the unsteady base-state in such a way that 
the perturbations and the base-state vary with time simultaneously (see Appendix [A] for the 


algorithm used). For this purpose, Eqs. (21)-(23) are solved using a highly accurate pseudo- 
spectral method to convert the linearized perturbed equations into a system of ordinary 
differential equations with algebraic constraints. We apply discrete Fourier transform to all 
the unknown variables, 


c'(x, y,t) = Zv,,(t)e <k ’ x+k ' y) , 


p,g 


a,m 


,y,t) = ^2cip,q{t)e 


i(k v x+k q y) 


P.9 




(26) 

(27) 

(28) 


P.9 


We also consider the discrete Fourier transform of the multiplicative terms, 

dc b d'll’ 1 


N(x,y,t) = 


dx dx 




,i(k p x+k q y) 


P.9 


j{x,y,t) = = Y, J P , q (ty (kpX+kqy \ 


dx dy 


P.9 


Mx,y, t) = = Y. t f .,,(«)e‘ ( ^"+‘.»>, 


dx dy 


(29) 

(30) 

(31) 


P.9 


where k p = 2 r Kp/L x ^ k q = 2Tcq/L y , p, q G N U {0}. The coefficients of the Fourier trans¬ 
forms (c Pi9 , C\ pcp etc.) are computed using the fast Fourier transform (FFT) when¬ 
ever c'(x,y,t), c' a m (x,y,t), i^'(x, y, t), etc. are known at the collocation points, x p = 
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ip/N x )L x , p = 0,1,..., N x - 1 and y q = ( q/N y )L y , q = 0,1,. .., N y - 1. Here N x and 
N v correspond to the number of spectral points in the longitudinal and transverse direc¬ 


tions, respectively. Substituting Eqs. (26)-(31) into Eqs. (21 )-(23) following differential 
algebraic equations are obtained, 


d t 

ddip q 

d t 


' + 5Jlp,q 


^p,q = R[Nj 




( k p + kg) &p,qi 

= - [i\k p + 5 (fc 2 + k 2 q )] ci Pi „ 
+ ik,c p ,„) / (kl + k\) . 


(32) 

(33) 

(34) 


Operator splitting method is employed to solve the differential equations (Eqs. (32) and 


(33)), subject to the algebraic constraints (Eq. (34)). Using the values at the time level t, 
we predict Cp i9 (i + At) and c\ pq (t + At) using the second order Adams-Bashforth method 
which are corrected using a trapezoidal rule. Inverse fast Fourier transform is used to obtain 
the corresponding values, d(x , y,t + At) and c' am (x , y,t + At), in the real space at the next 
time step t + At. Multiplicative terms are calculated in the real space at the new time level 
(the detailed algorithm of this numerical method can be found in Tan and Hornsy [8]). 

In order to use the benefit of the FFT we employ periodic boundary conditions on both 
the longitudinal and transverse boundaries, which are obtained straightforward from the 


physical boundary conditions, Eqs. (24) and (25). Since, x = 0, £ = 0 is a singular point 


for the error function base-state concentration profiles (Eqs. (15) and (16)), numerical 
simulations are performed by taking the initial time to = 1CU 3 . Convergence study has been 
carried out by taking spatial discretization steps (Ax, Ay) = (4,8), and (Ax, Ay) = (4,4) 
in a computational domain [0,4096] x [0,512], Relative error has been computed using the 
standard Euclidean norm on M 2 for amplification measure (see Eq. (38)) and it is found 


to be of 0( 10~ 2 ). To get optimal result, thus (Ax , Ay) =(4, 4), with At = 0.1 has been 
chosen. In the next section, the amplification measure and growth rate of the infinitesimal 
disturbances are presented. 
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FIG. 2: A straight line fitting to the normalized amplification: (a) In (E m y(t)/E m y(t 0 )) 
and (b) \n(E c v(t) / E c v(t 0 )) for k' = 0.2, R = 3 ,w = 512 ,k = 0.1. 


D. Amplification and growth rate of the perturbations 

We use sinusoidal perturbations of the form, 

{ e * cos (ky), at x = x r and x r + w 

(35) 

0, otherwise, 

where x r is the position of the rear interface of the finite slice, k is the non-dimensional wave 
number, e is the amplitude of the perturbation, which is taken as 10“ 3 and g corresponds to 
c', c' a m , and ?//• The present I VP approach measures the growth rate of each flow parameters 
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individually which are obtained from amplification measures, without invoking QSSA. To 
quantify the amplification gain at time t, it is necessary to define a norm which is generally 
based on the kinetic energy of perturbations. We define 

E cV {t) := [ (c'(x,y,t)) 2 + (u'(x,y,t)) 2 dxdy, (36) 

Jn 

E mV (t)\= [ (c' a m (x,y,t)) 2 + (u'(x,y,t)) 2 dxdy, (37) 

Jn 

where is the computational domain. E c y and E m v correspond to the amplification mea¬ 
sures associated with the solvent and solute concentration perturbations coupled with the 
velocity perturbations, respectively. Similarly we can define, E c (t) := f Q (c'(x, y, t)) 2 dA and 
Emif) : = fn( c a,m( x, y, t)) 2 dA to quantify the amplification measure of the solvent and solute 
concentration perturbations individually. 

The temporal evolution of the logarithm of the normalized amplification measures E c y(t) 
and E mV (t ) are shown in Fig. ([2]) for R = 3, w = 512, k = 0.1 and k' = 0.2. It is ob¬ 
served that both l\\{E c y(t)/E c v{to)) (see Fig. ga)) and \n(E mV (t)/E mV (t 0 )) (see Fig. gb)) 
increase linearly at early times, determining an exponential growth of the perturbations. 
Thus, following Kumar and Hornsy [26] the growth rate <7f(t) can be defined as follows. 


Definition III. 1 (Growth rate and growth function). The instantaneous growth rate of the 
perturbations is defined as, 

_ ld[ln(£/(t))] 

= 2-d^-’ (38) 

and the growth function is defined as the product of the instantaneous growth rate and the 
the corresponding time, i.e., af(t)t, where / G {c,m,cV,mV}. 



E. Direct numerical simulations 


In this section we discuss direct numerical simulations of the perturbation equations using 
the pseudo-spectral method [8j. Unlike LSA, in this case a finite amplitude perturbation is 


allowed by substituting 9 — 1 in Eqs. (18)-(20). The resultant equations are solved following 


the same algorithm discussed in Sec. Ill C The nonlinear (multiplicative) terms are defined 
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as, 


dc b d'lp' dd d't/j' dd d'd' 
dx dx dx dx dy dy 


Np,q{t)' 


t i(k p x+k q y ) 


( 39 ) 


P,Q 


dc b Oil)' dd dil>' dd dtl>' 
dx dy dx dy dy dx 


Jp,q(t) 


e i(k p x+k q y) 


(40) 


p,q 


T( K\,mdi)' dd d^,' dd dd' 

Ji{x,y,t) = — -— + 


dx dy dx dy 

J2 Ji P , q (t)d^ x+ ^ y) . 


dy dx 


(41) 


p,q 


For the validation of the numerical code the results of Mishra et al. |19] are reproduced. 


Further, the growth rate obtained from DNS is calculated from Eq. (38). In the absence of 


adsorption, i.e., k' = 0, the solute remains in the mobile phase and follows the dynamics of 
the sample solvent and acts like a passive scalar. Thus, the growth function <r mV t associated 
with the solute concentration must be the same as that of the solvent concentration a c yt 


see Eqs. (18) and (20)). 


IV. STABILITY ANALYSIS WITHOUT ADSORPTION 


A. Displacement of two semi-infinite fluids 


In this section, the linear stability analysis of classical VF instability between two miscible 
fluids in a porous medium is revisited in the absence of adsorption, i.e., when k' = 0. It is 
observed that letting w —> oo in Eqs. (15) we obtain c b (x,t ) = 1/2[1 + erf(a:/2V^)], which 
corresponds to the stability problem associated with the displacement of two semi-infinite 
fluids [ID] . In the framework of the present formulation this can be obtained by considering 


the sample width w very large in Eqs. (15) and (35), so that the right interface remains 


unaffected from the dynamics at the left interface. Corresponding stability equations are 
solved using the IVP method described in Sec. |mc| 


Following the definition in Eq. (38), growth functions, a c t,ayt and a c yt are calculated, 


where ay is measured from Ey = E c y — E c . The temporal evolution of these growth functions 
are shown in Fig. [3] for R = 2, k = 0.1. This figure depicts that all the disturbances initially 
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FIG. 3: Temporal evolution of the growth functions for two semi-infinite fluids with 

R = 2,k = 0.1. 


decay in diffusion dominated regime and grow when convection dominates, and eventually 
become unstable at a later time. It is observed that for all t > t 0 , a c t < a c yt < ayt, 
which is consistent with the results of Tan and Homsy |T0|, who discussed only a c and ay 
in their study. Further, it is verified that perturbations to each flow variables are essential 
in hydrodynamic stability analysis. We observe that if perturbation is not imposed to ip', 
ay becomes unconditionally unstable for all t > t 0 , which violates the physics of early time 
diffusion dominated regime. As the velocity plays an important role in VF instability, the 
growth function a c yt is used to analyze the onset of instability for the classical VF problem 


and also for finite slice displacement (see Sec. IVB). The onset of instability is defined as, 


t c = min {r > to : a c y(t < r) < 0, a c y(t > r) > 0} . (42) 

T 


In summary, the present LSA method differs from the other existing linear stability 
analyses by two conceptual distinctions. It successfully captures both the onset of instability 
and the early time diffusion dominated regime, which was never achieved with the well known 
QSSA methods [10i, 12] fl3] . Also, the present LSA measures the magnitude of the growth 
function is close to zero at the early time, which represents that at such initial period the 
growth function is of order 0(1) [ID]. 
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(a) 



(b) 



FIG. 4: (a) Temporal evolution of the growth function, (J c yt. for R = 2, k = 0.2 and 
different width of the finite slice, w. Symbols correspond to growth function associated 
with w = 128(0); 256(D) and oo(0). (b) Magnified of (a) near the onset of instability. 

B. Finite slice displacement 

Here we are interested to investigate the influence of the finite slice width w on the 
onset of VF dynamics between the displacing fluid and the sample solvent. Fig. [4] depicts 
the temporal evolution of a c yt for R = 2, k = 0.2 and different sample width w. For these 
parameters the critical sample width for the onset of instability has been found to be w c ~ 10 
(see Fig. [l]). It is shown that the onset of instability are indistinguishable for all w > w c 
(see Fig. |4](b)) . Further, we observe that the growth functions corresponding to a finite 
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slice of width w > 64 (see Fig. |4](a)) and the displacement of two semi-infinite fluids are 
indistinct. This result is consistent with the linear stability results of Pramanik and Mishra 
[13] and numerical simulations of De Wit et al. [[27]. For the increasing values of R the 
viscosity contrast increases which influences the disturbances to grow faster. Therefore, it 
can be shown that the critical sample width w c decreases as R increases m- 



FIG. 5: Growth functions, a c yt ( - ) and a m vt ( 0 )> f° r 

k' = 0, R — 2, w — 512, k = 0.15. The solid dot (#) represents the onset of the nonlinearity. 

Inset: Magnified near the onset of instability. 


In the absence of adsorption of the solute on the porous matrix, the present model reduces 
to the model of De Wit et al. [27], and this was confirmed by Mishra et al. [T2] through 
DNS. Before we analyse the influence of retention parameter, k', on the stability of the 
adsorbed solute, we confirm the same through the present LSA. In this context we choose 
R = 2, k = 0.15, k! = 0, and compare the temporal evolution of the growth functions, a c yt 
and cr m yt. These growth functions are shown in Fig. [5j which depicts that a c yt and a m vt 
are indistinguishable, thus confirming the validity of the present IVP based LSA. Further, 
the growth rate obtained from DNS is compared with those obtained from LSA in Fig. [5} 
It shows that the DNS results coincide with the LSA results before non-linearity becomes 
dominant at t ~ 167.5. Therefore, linear stability theory is not valid beyond this time and 
one needs to solve the complete nonlinear problem to capture the dynamics of the instability 
pattern. 
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V. STABILITY ANALYSIS OF AN ADSORBED SOLUTE 


In this section the numerical results obtained from both LSA and DNS are discussed to 
analyze the effect of the retention parameter, A, on the growth of the perturbations associ¬ 
ated with the solute dynamics. As k' has no effect on the temporal evolution of the sample 
solvent c, we choose the growth function cr m y£ associated with the solute concentration c Qim 
to quantify the effect of linear adsorption isotherm on solute dynamics. 



k 

FIG. 6: Neutral curves for R — 2, w — 512 and different values of the retention parameter 
k' = 0, 3 and 8. The critical time t c and corresponding wave number k c on the parameter 

space is marked by solid dot (#). 


A. Neutral curves and dispersion curves 

The overall characteristic of the temporal instability driven by viscous forces can be 
compactly presented in a phase space spanned by the perturbation wave number k and time 
t. Fig. @ represents the isocontours of the growth rate, cr m y = 0, in the (k,t) parameter 
plane for R = 2, w = 512 and the retention parameter k' — 0, 3 and 8. The area above each 
contour line corresponding to different k! represent the unstable region for the respective 
values of the retention parameter. For small times all perturbation wave numbers are stable, 
later a band of wave numbers become unstable. The critical point ( k c , t c ) is marked with 
a solid dot at the lowest point on the isocontours a m y = 0. We show that the region of 
instability is larger for non-adsorbed case (V = 0) in comparison to the case when k' ^ 0. 
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This signifies that adsorption of the solute concentration on the porous matrix reduces the 
instability. Further, it is observed that the onset time t c has a non-monotonic dependence 
on k'. The detail analysis of this non-monotonicity and the underlying physics are discussed 
in Sec. EJ3 



FIG. 7: Dispersion curves at different times for w = 512, R = 2: k' — 0(0), 3(D), 8(0)- 


Dispersion relation is a convenient measure of the growth rate corresponding to the differ¬ 
ent wave numbers of the disturbances. Fig. [7] depicts the dispersion curves corresponding to 
the same parameters as in Fig. [6] at different times. This figure shows that adsorption of the 
solute on the porous matrix has an overall stabilizing influence. Physically, this illustrates 
the fact that fingers which are too narrow are immediately smoothed out by the increase in 
transverse dispersion. Not only the growth rate of the perturbation associated with certain 
wave number is reduced with increasing k 1 , the spectrum of the unstable wave numbers is 
also reduced as the retention parameter is changed from k' = 0 to k' ^ 0. Besides showing 
the spectrum of unstable wave numbers, dispersion curves also represents the magnitude of 
the growth rate of the unstable modes. At a given time the maximum possible instanta¬ 
neous growth rate that can be achieved by any initial condition is denoted by cr™y-(t) and 
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is defined as, 


<CvW= SU P <Tmv{k,t), (43) 

0<k<oo 

and the corresponding wave number is the most unstable wave number, which is denoted by 
fc max (t), at time t. Fig. [7] shows that for R = 2,w = 512 and t < 10 the most unstable wave 
number, fc max , which decreases with t, is almost independent of the retention parameter, k' 
(see the corresponding curves for k! = 0,3,8). Corresponding cr™y x increases with t having 
a non-monotonic dependence on the retention parameter k', for t < 6. 



FIG. 8: Dominant wave number, /c max , of the solute concentration for w = 512, if = 2 and 

different values of k'. 

In Fig. [8] the temporal evolution of the most unstable wave numbers are presented for 
R — 2,w — 512 and the retention parameter, k' = 0,0.5,1 and 2. This figure depicts 
that for all k fc max decreases monotonically with t, indicating stabilization of short wave 
perturbations at later times. We show that, fc max is almost independent of k' at early times, 
while at later times k max decreases monotonically when k! increases, which signifies that the 
long wave perturbations are smoothed out more rapidly with increasing k 1 . 


B. Optimal growth 

In an experiment or a real physical system a perturbation consists of combination of 
different wave numbers, so it is important to calculate the onset of instability considering 
all possible wave numbers. Fig. [9] represents for R, — 2,w — 512 and k! = 0,3,8. 
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FIG. 9: Optimal growth function, cr™yt, of the solute concentration for w = 512, R = 2 
and different values of k' = 0, 3, 8. Inset shows the non-monotonicity of onset time and 

dominance of diffusion at early times. 


It confirms that adsorption of the solute on the porous matrix results into an overall delay 
of the onset of instability. The qualitative as well as quantitative effect of k' on the onset 
of instability are discussed below. It must be noted that the retention parameter has no 
effect on the dynamics of the sample solvent. From the mass-balance equation for solute 
(see Eq.Q) it is seen that the dispersion coefficient is 1/(1 + n!\ while the travelling wave 
velocity of the mobile phase solute concentration is — n '/(1 + fd) in the longitudinal direction. 
Therefore, as k' increases the solute moves in the upstream direction away from the sample 
solvent more rapidly and eventually disengages from the solvent zone. In this process the 
rear or the frontal interface of the solute zone features VF instability depending on the value 
of the retention parameter, k' [12j. As n! increases, from Eq.Q it can be observed that the 
dispersion coefficient decreases and hence the optimal growth of disturbances for n' = 0 
damped most in comparison to n' ^ 0. This is depicted in the inset of Fig. [9j Further, for 
n' > 3, the onset time is early in comparison to when k' G (0, 3]. Thus the onset time varies 
non-monotonically with respect to hi'. 
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FIG. 10: Onset time, t c calculated from a m v, as a function of the retention parameter, k 1 , 

for R = 2,3,w = 512, k = 0.1. 


C. Effect of retention parameter and log-mobility ratio on the onset of instability 


In order to quantify the influence of the retention parameter k' on the onset of instability, 


we solve the IVP, Eqs. (21)-(23), for two different log-mobility ratio R = 2 and 3, sample 


width w = 512 and k' G (0, oo). Thus calculated onset time t c is shown as a function of k' 


in Fig. 10 This figure represents the onset time calculated from a mV and depicts that t c 


depends non-monotonically on k'. In order to analyze the non-monotonic behavior of the 
onset of instability more precisely, for a given wave number k and log-mobility ratio R we 
define the largest onset time R ), k, R) over all retention parameter k' as, 


£“(«!, k, R) = max t c (k, R), 


(44) 


and the corresponding retention parameter is denoted by k) ( R). In particular, we choose 
k — 0.1 to compute t r ( and the corresponding n J(f?) for two values of log-mobility ratio R — 2 
and 3. It is determined that t^{R = 2) = 6.26, k[(R — 2) = 3 and tf(R = 3) = 2.81, k^R = 


3) = 1 and they are represented by solid dots in Fig. 10 More generally, it can be shown that 


.Ri < f? 2 implies k[(Ri) > k , 1 (R 2 ) and tf(K[(Ri), R±) — t c (0,i?i) > R 2 ) ~ t c (0, R 2 ). 

We further show that, for sufficiently large values of k', the instability in the adsorbed 
solute sets in earlier than the solvent. The corresponding value of k’ depends on the log- 
mobility ratio, and it decreases with decreasing R. For large values of k! the axial dispersion 
of the solute concentration is smaller than that of the solvent concentration. Thus the 
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stabilization effect of the dispersion is more on the perturbations of the solvent concentration 
in comparison to the solute concentration. Hence, in the early time diffusion dominated 
regime a m decays slower compared to a c . As a consequence, <j m y becomes larger than cr cV 
and hence shows an early onset for the solute concentration than the solvent concentration. 



t 

(b) 



FIG. 11: For R = 2,w = 512: (a) Contribution of fingering to the standard deviation cr a j 
as a function of time, (b) Onset time, if, obtained from DNS as a function of the retention 
parameter, hi'. Inset: Density plot of the solute concentration field for n' = 0.2 and hi' = 5 
at the corresponding onset time, if = 740 and t f = 640, respectively. 

Next we analyze the dependence of the onset of fingers in the nonlinear regime on n! by 
solving the fully nonlinear equations [19:, 27]. De Wit et al. [37] showed that the standard 
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deviation of the solvent concentration zone due to viscous fingering <7/ starts growing from 
zero at the onset of VF. Further, following their analysis, Mishra et al. j!9| measured 


</(*)= rL 


r^ x cIt 

JO b a,m 

f Lx c dx~ \ r 

Jo Ca ’ m \ Jo 


rL_ 

Jo ~ XC ' 


a,m 


dx 


Ca,m dx 


w 2 2 1 

~ 12 + 1 + K' 


(45) 


for the solute concentration zone, where c a ,m{x,t) = f^ y c atTn (x, y, t)dy is the transversely 
averaged solute concentration profile. For R = 2,w = 512, Fig. [TT^a) depicts the contri¬ 
bution to the standard deviation due to viscous fingering, d a j , as a function of time for 
k' = 0, 0.5, 3. From this figure, it is seen that the time, when d a j starts deviates from zero, 
varies non-monotonically with respect to k'. Such a time was mentioned as the onset of VF 
by Mishra et al. [IHJ. In order to quantify this precisely, we define the onset time of fingers 


as, 


= min {t > t 0 : d a j(r > t) = 0, cr a j(r < t) > 0} . (46) 

T 

In Fig. [TT|(a) solid dots (#) correspond to the onset time measured by visually observing fin¬ 
gers from the density plots of the solute concentration. For R = 2, w = 512 the dependence 
of as a function of k' has been plotted in Fig. 0 b ). which represents a non-monotonic 
dependence of tf on k'. It also shows the density plot of solute dynamics (see inset of Fig. 
[TTfb)) at the onset of fingers, It can be shown that this non-monotonic characteristic of 
onset of instability can not be captured if t c is calculated from a m , instead of cr m v- Thus, 
we conclude that the growth rate measured from the amplification measure of the solute 
concentration coupled with the velocity perturbation is more realistic than that correspond¬ 
ing to only solute concentration. Further, the effect of the log-mobility ratio R is consistent 
with the classical VF instabilities in two component models mm, e.g., for a fixed n' ^ 0 


and wave number, instability sets in earlier for larger R (see Fig. 10). 

Mishra et al. [19j calculated that in absence of viscosity contrast between the displacing 
fluid and the sample solvent (i.e., R — 0), the disengagement time of the solute from the 
solvent is given by (see Eq. (22) in [HI]) 

tcric = (1 + 1/AC ) 2 \/2(l + \J8/(1 d~ V)) + J 2(1 + \/5/(\ + V)^) + k! wj (1 + V) . (47) 


2 

In the limiting case of very large V, t cr i c attains an asymptotic value (\/2 + \J 2 + w)~. As 
an example, for a finite slice of width w = 14 this asymptotic value becomes ~ 29, which is 
very large compared to the onset of fingering instability t c measured from the present LSA. 
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(a) 



(b) 



FIG. 12: (a) Growth functions, <v m yf, of the solute concentration for 
R = 2, w — 512, k — 0.1 and k' = 0 (solid line), k' = 0.2 (dashed line): Lines without 
symbol correspond to LSA and lines with symbols represent DNS. (b) Magnified of (a) 
near the onset of instability marked by the solid dots (#). 


Thus the disengagement of the solute from the solvent zone is not possible before the linearly 
unstable modes set in for instability. However, onr LSA captures the non-monotonicity in 
t c . Thus we conclude that for all k' ^ 0 there exist linearly unstable modes, which may not 
develop into fingers in the nonlinear regime for large k'. 

Next, we investigate the influence of k' both in the linear and nonlinear regimes and 
compare the results obtained with those in the absence of adsorption. In this context we 
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choose R — 2,w — 512, k = 0.1 and k' = 0,0.2. The temporal evolution of the growth 


functions obtained from both LSA and DNS are shown in Fig. 12 Fig. 12 b) illustrates 
that when diffusion dominates at the early time, the growth functions cr m yt corresponding 
to two values of k' are almost identical (see curves for t < 2). As soon as they start growing 
in the convection dominated regime they are different (see curves for t > 2), and this leads 
to different onset time. However, it is identified that the onset of non-linearity is almost 
independent of the retention parameter k' (see Fig. |l2)(a)). 


VI. CONCLUSION 

We have theoretically studied the onset of fingering instability in a finite slice of linearly 
adsorbed solute. Instability is driven by the viscosity contrast between the displacing fluid 
and the sample solvent containing the solute. We presented a linear stability analysis based 
on a Fourier pseudo-spectral method, which, compared to QSSA methods, captures both the 
onset of instability and the early time diffusion dominated regime. The linearized equations 
are solved as an initial value problem and the growth rate associated with each perturbation 
quantities is calculated from their respective amplification measure. The numerical results 
revealed that the exponential growth of the perturbations are reasonable in the linear regime. 
It is verified that, in absence of retention, i.e., k' = 0, the onset time is independent of 
w > w c , a critical value. Further, it is shown that there exists a threshold finite slice 
width, beyond which the stability analysis for finite sample is identical to that of a single 
interface displacement. Another very interesting observation is that, the onset time is a non¬ 
monotonic function of the retention parameter, k'. It is shown that for a given wave number 
the largest onset time and the associated retention parameter decrease as log-mobility ratio 
increases. The present LSA agrees qualitatively with DNS and we successfully distinguish 
between the linear and nonlinear regimes. Analysis with velocity dependent dispersion and 
nonlinear adsorption isotherm has been undertaken for further study. 


ACKNOWLEDGEMENTS 

S.P. acknowledges the National Board for Higher Mathematics, Department of Atomic 
Energy, Government of India for the Pli.D. fellowship. 


25 




Appendix A: Algorithm of the present IVP approach for LSA and DNS 

The present linear stability analysis is of generic type as it handles the unsteady base- 
state very carefully that helps to capture the underlying physics more appropriately. Below 
we describe the algorithm of the numerical method used in the present LSA: 


1. We introduce perturbation to c, c a ,m and -0 at time t 


in Eqs. (21 )-(23) at t — to- 


to and compute the coefficients 


2. Time integration is performed by taking the perturbations as the initial condition to 
the unknown variables, c', d am and ^. 

3. Obtained solutions, c',c' arn and ip' are used as the initial condition for the time march¬ 
ing in the next step. Repeat this step until desired result is obtained. 


Appendix B 


Most of the LSA methods existing in the literature of fingering instabilities driven by 
viscosity contrast solve the following linearized equations, 


3 2 d 2 dc b d ' 
U-- 


dx 2 dy 2 dx dx 

8 _ d 2 _ d 2 
dt dx 2 dy 2 


v! = -R 


d 2 c’ 
2 ’ 


dy 


c = 


dc b , 
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dx 


—-h—-h— A — 
dt dx 2 dy 2 dx 


dc h 


/ r ~ a.m / 

c„ m = — 0 ——u . 


a,m 


dx 


(Bl) 

(B2) 

(B3) 


Since, the coefficients of Eqs. (Bl)-(B3) are independent of y, wave like disturbances are 
assumed of the form, 


(“'.c'AmlO'ii/.t) = (4 , , < *,0)OM)exp(!%), 


(B4) 


where k represents the non-dimensional wave number in y-direction. The resultant linear 
equations can be written compactly as a non-autonomous IVP 


dq (t) 

dt 


= A(t)q(t), q(x,t 0 ) = q 0 (x), 


(B5) 
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where 
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f)r b dd dc b 

M 4 = - —, M 5 = M? - A-^, M 6 = 

ax ax ax 


I is the identity operator, and t 0 is the time at which the perturbations are introduced. 


1. Quasi-steady state approximation method 


This classical frozen time approach of investigating the instability assumes that the un¬ 
steady base state evolves very slowly in comparison to the perturbations. Rewriting Eq. 


(B4) as (u\ d, d )(x, y, t) = (T, <h, @)(x) exp (iky + a(t q )t), where t q is the time at which 


the unsteady base states (see Eqs. (15) and (16)) are frozen, Eq. (B5) reduces to an algebraic 
eigenvalue problem, i.e., A(t) q(x) = cr(i ? )q(x). The maximum eigenvalue of (A-\- A t )/2 can 
be interpreted as the maximum possible instantaneous growth rate [28] that can be achieved 
by any initial condition at early times. This reveals that QSSA does not capture separate 
growth rates for the solvent and solute concentration perturbations, thus restricting the 
analysis of present model to that of classical VF instabilities [T21 133] • However, the DNS 
results of Mishra et al. [T9] showed that the retention parameter k’ influences the instability 
of the solute, but not the solvent. Therefore, the underlying instability dynamics of solute 
can not be captured by the QSSA method. Hence, an IVP approach is imperative to find 
the LSA of adsorbed solute transport. 


2. IVP approaches 


Due to unsteady nature of the base state flow the growth rate obtained by solving the 


IVP, Eq. (B5), is sensitive to the initial condition q(x,£o) = Qo^)- I 11 the literature, the 
IVPs are solved using an initial condition that corresponds to a random perturbation in the 
whole spatial domain, 


q(x, to) — e * rand(x), Vx, 


(B6) 
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where rand(a:) represents a random number generator between —1 and 1, and e corresponds 
to the amplitude of the perturbation, which is a very small positive number (0 < e 1). 
Since it is known that the fastest growing perturbation is localized around the diffusive 


interface [9j, we solve Eq. (B5) in ( x,y,t ) coordinate system with an initial condition of the 
form, 


qOMo) 


{ e * rand(x), a\ < x < bi 
0, otherwise. 


(B7) 


Here the spatial interval [ai,&i] corresponds to the thickness of the diffusive layer, and its 
location depends on the viscosity contrast between the fluids. The results obtained have 
very good agreement with the DNS results and the present pseudo-spectral method based 
LSA. 
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